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Abstract 

We insert some asymmetries in the continuous Abelian sandpile models, such as directedness 
and ellipticity. We analyze probability distribution of different heights and also find the field 
theory corresponding to the models. Also we find the fields associated with some height variables. 
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1 Introduction 

After the seminal work of Bak, Tang and Wiesenfeld [I], sandpile models have been found to 
be a very helpful framework to study Self-Organized Criticality (SOC). SOC is believed to be 
the underlying reason for the scaling laws seen in a number of natural phenomena [2]. Further 
investigations showed that BTW sandpile model has an Abelian property and was renamed Abelian 
Sandpile Model (ASM) [4j. The model is still the simplest, most studied model of SOC, in which 
many analytical results has been derived. For a good review see ref. [3]. 

Extensive work has been done on this model, analytical and computational. To name a few, one 
can list the works by Dhar [HE], where the probabilities of some height cluster were computed. In 
[3], Dhar computed the number of recurrent configurations and showed that all occur with equal 
probability. Later, in [5], Majumdar and Dhar calculated the probabilities of occurrence of some 
specific clusters, known as Weakly Allowed Clusters (WAC's). The simplest of these clusters is one- 
site height one cluster. The probabilities of other one-site clusters with height above 1 was computed 
in ref. [6]. There are many other analytical results, among them one can mention the results on 
boundary correlations of height variables and effect of boundary conditions EE EH EE] 12] 5 on 
presence of dissipation in the model [9] [TO], [J3] , on field theoretical approaches [9] [HI [JS] [16] [T7] , 
on finite size corrections [S][IT] and many other results [3j. 

In ASM, the height variables take only integer values. However there are some other models 
in which the height variables can take real values [181 [T9~] [2"T| [22] [23] . In Zhang model, which is 
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not Abelian, when a site topples all of its energy (height) is distributed among its neighbors. On 
the other hand the model introduced by Gabrielov, called Abelian Avalanche Model (AA Model) 
and the model introduced by Ghaffari et. al. have the Abelian property. In [23] a continuous 
version of ASM is considered on a square lattice and many of its different properties are studied. 
When we let the heights be real number, naturally the components of toppling matrix could be real 
too. In this paper we consider continuous ASM [23J and change different components of toppling 
matrix. We will do it in a way to add spatial asymmetries to the model. Two specific asymmetries 
are considered, one is to introduce a preferred direction. The second is to introduce elliptical 
asymmetry. The first one is related to directed sandpile model and the second, as we shall see 
results naturally from the original BTW model. We have derived the probability distribution of 
height variables and also avalanche distributions. Also the field theory associated with these models 
are considered. The new models are still critical and show scaling relations. 

In the next section we first review the Continuous ASM. Next we consider the effects of a 
perturbation that introduces a preferred direction in the model. In the forth section we introduce 
ellipticity in the model and derive many properties of such a system. 

2 Continuous Abelian Sandpile Model 

The Continuous Abelian Sandpile Model (CASM) on a square lattice is defined in [23]. Consider a 
L x L square lattice. To each site u = a continuous height variable, h(u) is assigned, where 
without loss of generality we assume these variables are in [0,4) (in the original paper it was taken 
to be in [0,1) ). The evolution rules consist of the following rules: 

1) At each time step, a site is selected randomly and an amount of sand is added to it and other 
height variables are unchanged. The amount of sand added to the site is a random real number in 
the interval [p, q] £ [0,4). The distribution function could be any function, but for simplicity we 
take it to be uniform on the interval [p,q\. If height remains below four, the new configuration is 
stable and we go to the next time step. 

2) If height of the sand at that site becomes equal to or greater than four, the site becomes 
unstable and topples in the following way: it gives an amount of sand with height one to either 
of its neighbors. In other words h(v) — ► h{u) — A c (u,v) for all v, where A c (u, v) is the toppling 
matrix of the CASM and is defined as 



So if the height of a site is 4 + a, after toppling its height will be a and the heights of neighbors 
will be added by the amount 1. As a result of this toppling, some of the neighbors may become 
unstable and an avalanche may occur. 

It has been shown in [23J that there is a (many to one) mapping from configurations in CASM 
to configurations in ASM, which preserves the dynamics, therefore the new model reproduce most 
of the interesting features of ASM. Also it is shown that the probability distribution is piecewise 
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Figure 1: The probability density profile of height variables 



constant. The probability of finding a site with height h £ [ft — 1, ft) with ft = 1, 2, 3, 4 is equal with 
p(ft), the probability of finding a site with height ft in the usual ASM. The claims were proved and 
were in agreement with simulations. 

As the heights are continuous, we can take other integer parameters of the model to be real. 
This has been done to the dissipation parameter in the model[23j. The probability distribution is 
found when there exist dissipation with different (real) values in |23j . Yet there are some other 
parameters which are integer in the original model, but may be taken to be real in the new model. 
In this paper we consider the amount of sand transferred to the neighbors and let it take different 
real values and see how the properties of the model are changed. 



3 CASM with preferred direction 

One of the easiest modifications to the toppling matrix is to introduce a preferred direction, that is; 
sand grains are more likely to move to left rather than right. Such modification is also studied in the 
usual ASM [23], but as the toppling matrix components have integer values, the only possibility is 
to send two sand grains to left and nothing to right. The same can be done in the vertical direction. 
If this is the case, then the toppling matrix will be upper triangular and all the configurations will 
be recurrent [22]. However, in our model we are able to control the amount of directedness: when 
a site is toppled, then 1 + e grains of sand move to left and 1 — e grains of sand move to right. 
Clearly e = corresponds to the usual ASM and e = 1 corresponds to the case mentioned above, 
provided we add both horizontal and vertical directedness. 

Several properties of the model can be studied, the probability distribution, the corresponding 
action, the avalanche distributions etc. The probability distribution of finding different heights is 
sketched in Figure 1 for some different values of e when both directions are directed. As we see, 
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still we have the same step-like pattern as in the case of undirected CASM. Some new small steps 
have appeared which have the length proportional to e. As e becomes larger, these steps become 
longer and little by little the original pattern fades away. At e = 1 probability distribution becomes 
uniform, as it is expected [23] . 

Also the avalanche distribution is studied in this model. System sizes from 64 to 1024 has been 
investigated. Starting with a lattice of randomly distributed heights h £ [1,4), the system is evolved 
to reach steady state. After that, we begin the measurements. The measurements are averaged 
over about 10 7 avalanches. We studied two characteristics of an avalanche (though we bring the 
result of one of them) : the total number of toppling events is called the size s of an avalanche and 
the number of distinct toppled lattice sites which is denoted by Sd, ■ Because a particular lattice 
site may topple several times, the number of toppling events exceeds the number of distinct toppled 
lattice sites, i.e., s > Sd, except when the system is fully directed in both directions in which case 
we have s = Sd- In the critical steady state the corresponding probability distributions should obey 
power-law behavior characterized by exponents t s and r^: 

P(s)ocs- Ts , (2) 
P(s d ) oc s d T « (3) 

Figure 2 displays the obtained results for the distribution P(s) for different system sizes and 
three different values of e = 0.1, 0.4, 1.0 when the both directions are directed. A powerdaw fit 
to the straight portion of these curves yields the exponents t s (L). Figure 3 shows a plot of the 
exponents r s (L) vs 1/logL for e = 0.1 and 0.4. This allows us to extract the exponent r s (oo). As 
example, for the cases shown in Figure 3, it turns out to be: r s (oo) = 1.31 =t 0.08 for e = 0.1 and 
t s (oo) = 1.38 db 0.08 for e = 0.4, though the error bars are quite big and we do not have reliable 
results at hand. It is worth to mention that for small values of s, the magnitude of slope is less 
comparing with the one associated with greater values of s which is due to finite size of lattice. It 
is seen that still the system is critical and self-organized. There is no characteristic length in the 
model and finite size scaling is clearly observed from Figure 2. The exponent r s seems to increase 
as e becomes greater which could be expected from the known behavior of directed sandpile model. 

In the continuum limit, we may assign an action to the theory. As in the steady state, only 
recurrent configurations appear, and they appear with equal probability, the partition function is 
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Figure 2: Avalanche size distribution for e = 0.1, 0.4, 1.0 and for lattice sizes L = 
64, 128, 256, 512, 1024. 
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Figure 3: The exponents t s (L) is a linear function of 1/logL. The intersection with vertical axis 
gives t s (oo) 

simply equal with the total number of recurrent configurations, that is, the determinant of toppling 
matrix. Such determinants could be written in terms of integration over grassman variables: 



In this way the action of the theory may be obtained. If lattice spacing is small, one can go to 
continuum limit and find a field theoretic action for the theory. In the case of usual ASM the action 
turns out to be the c = — 2 logarithmic conformal theory action |25|, I14j: 



For the theory defined above, the action is more or less the same, but there are some added terms 
due to directedness: 



where a is a constant proportional to e with dimension of length. This new term grows under 
renormalization group, therefore the large scale properties of the model is given by the fully directed 
model. 

4 Elliptical asymmetry 

In the previous section, it was supposed that the lattice has a preferred direction. The field 
associated with this deformation was a relevant one and hence we expected the theory show different 
characteristics on large scales. Also the toppling matrix was not symmetric: the amount of sand 
transferred from site i to the site j is not necessarily equal with the amount of sand transferred 
from j to i. Therefore one can not apply the burning test to the model. However it is possible to 
modify the toppling matrix in the following way so that it remains symmetric and, as we shall see, 
the added term to the action would not be relevant. 

We assume that the vertical links can carry an amount of sand equal to 1 — e and the horizontal 
links can carry an amount of sand equal with 1 + e (0 < e < 1). In this way the total amount 
of sand removed from a toppled site is 4 and the mass of sand is conserved. As we see the new 
toppling matrix is symmetric, but in real space, x direction is quite different from y direction, in 
fact; we have added some ellipticity to the system and therefore it does not have full rotational 
symmetry, however this deformation preserves the scale symmetry. 
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Figure 4: The probability density profile of height variables in elliptical asymmetry 

If we go back to the original BTW paper [lj, we observe that elliptical asymmetry naturally 
could arise in their model. In their paper, BTW first define a one dimensional model and using a 
specific transformation, express it in some other parameters. The same is done to a two dimensional 
version of the model to arrive at BTW model, however it could easily be observed that if we apply 
the same transformation to the two dimensional model, we will not arrive at the BTW model, 
instead we arrive at a model with elliptical asymmetry. Let's explain it in more details. In the 
one dimensional model, an ordered array of heights Zi is considered. Sand grains enter from left 
and leave from right. If the difference of the heights of two neighbors become more than two a 
sand grain moves to right. The slope parameters hi are defined as hi = z% — Zi+i- Then as result 
of "toppling" at site i we will have hi — > hi — 2 and hi±i — > hi±\ + 1. The generalization to two 
dimensions would be to define the slope parameters as hi = 2zi — Zi t R — Zi t D, where z^r and z^r> 
are the height variables right and down to the site i respectively. The dynamics would be if the 
slope parameter of a cite (mean value of the slope there) is more than a threshold, one sand grain 
moves to right and one sand grain moves downward. The dynamics in terms of h variables would be 
hi -> hi- 4, hi,L;R;D;U -> hi^R;D;U + 2 and hi t RU;LD -» hi t RU;LD - 1, where hi t Ru-,LD are the slope 
parameters of right-up and left-down next nearest neighbors of the site i. As we see there is an 
asymmetry in the two diagonal directions and the system does not have full rotational symmetry. 

The probability distribution of finding different heights is sketched in Figure 4 for e E {0.1, 0.4, 1.0}. 
The simulation is done on a 1024 x 1024 lattice averaging over 10 7 samples. We still have major 
and minor steps when e is small. As e becomes greater, the minor steps become larger and at e = 1 
we have a two-step graph, which indicates that the model is now effectively a one dimensional one. 

Also the avalanche distribution is studied in this model. Just as the case of directed CASM, 
system sizes from 64 to 1024 has been investigated. Again, starting with a lattice of randomly 
distributed heights h E [1,4], the system is evolved to reach steady state. After that, the mea- 
surements are begun. The measurements are averaged over about 10 avalanches. Probability 
distribution of avalanche size is sketched in Figure 5 for e = 0.1, 0.4. The case of e = 1.0 is totally 
different, as it corresponds to one dimensional ASM. 

Using The same tools, the slope r s (oo) for different values of e could be obtained, for example we 
have T s (oo) e =o.i = 1-25 ± 0.06 and t s (oo) £= o.4 = 1-24 ± 0.06. Again system shows finite size scaling 
and is critical, but there is not a remarkable variations in the exponents like r s as e is changed. 
Possibly they belong to identical universality class. 
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Figure 5: The avalanche size distribution in a system with elliptical asymmetry. 
4.1 One-site height probability 

In CASM with no asymmetries, one is able to find the probability of having any height in the interval 
[0, 1]. This is done using the burning algorithm[5j: consider the recurrent configurations with the 
property that the height of the site i is in [0, 1]. As this site burns after all of its neighbors, it should 
be a leaf in the corresponding spanning tree. Therefore, the total number of such configurations 
could be obtained by counting the number of spanning trees in which site i is a leaf. In the modified 
model the same arguments holds for h £ [0, 1 — e), because the burning test still works and any site 
with h S [0, 1 — e) has to be a leaf in the corresponding spanning tree. Note that as we have links 
with weight 1 + e, the sites with h € [1 — e, 1 + e] may not be a leaf. 

Similar to the way we find the probability of one in ASM, we modify the the weights of the 
toppling matrix to ensure that the site i is a leaf. This modification of weights is only done on a 
finite number of lattice bonds, therefore the probabilities corresponding to these local restrictions 
can be calculated in terms of finite dimensional determinant. 

We modify the original model by removing the bonds to three neighbors of site i such that 
only one bond connects site i to the system. All four possible choices to cut the bonds are not 
equivalent, there are two different modifications: A\ and A2 shown in Figure 6. In the case A\, we 
delete all neighbor bonds of site i except a left or a right one which has weight 1 + e, and change 
the weight of this remaining bond to 1 — e to ensure the height of site i won't be more than 1 — e. 
In the case A2, one of the vertical bonds with the weight 1 — e remains unchanged and connects 
site i to system. 

The two new toppling matrices are of the form A^ 1 ) = A — B^ and A^ 2 ) = A — respectively 
where B^ and B^ are defect matrices: 

t t 



i i 

Al A2 
Figure 6: Two types of bond modification 
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(8) 



As we have left-right and up-down symmetries, the number configurations where hi 6 [0, 1 — e) is 
proportional to N = 2(det A^ +det A^ 2 )). So the probability of finding at most 1 — e amount sand 
at a site is given by: 

' ldet{I-B^G) +det(I - B^G), 



^"^detA 2 
where G = A -1 , is Green function matrix with the following integral form: 



Gij — G{ 



oo 



w dpi f w dp 



cos(ipi)cos(jpj) 



1 



2tt j — 2tt 4 — 2(1 + e) cospj — 2(1 — e) cospj ' 



(9) 



(10) 



The new toppling rules with elliptic asymmetry in CASM make the values of Green function to 
depend on the orientation of two sites an addition to their distances. Hence their calculation is a bit 
more tricky. For example, in the symmetric case it is enough to calculate Ga, all other components 
of G could be obtained easily thereafter. We have computed the Green functions analytically and 
a few of them which were necessary to derive P(l — e) are listed in Appendix. Using these Green 
functions, the determinants can be computed and P(l — e) is derived. The formula is very long 
expression and we do not bring it here, we have sketched it in Figure 7. The result is very good 
in agreement with the simulations done. As we see, for e = we arrive at the well-known result 
P(l) in BTW model. The probability gradually comes down till it vanishes at e = 1 which could 
be expected. It is also possible to calculate other one site probabilities such as the probability 
fe£ [1 — e, 1 + e). In these cases, you have some non-local constraints, just like higher-than-one one 
site probabilities in usual ASM. 

P(l - e) 




Figure 7: The probability P(l — e) as a function of e 
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Figure 8: The amount of rotational asymmetry of two point functions for e = 0.1 and 0.2. 
4.2 Two-point correlation function 

Again, using the same scheme introduced by [5], one is able to find correlation functions of finding 
sites with height less that 1 — e. Making the modifications according to configurations A\ and A2 
in the vicinity of any two sites, it is possible to find the probability that both sites have heights 
less than 1 — e. In order to compute this joint probability, we consider one site at the origin and 
another at i = (r, 6). There are four different cases depending on how we modify the toppling rules: 
the defect matrices associated with either of the site could be or B^ 2 \ Collecting all the four 
possible configurations, one finds: 



P(l-e,l- e ) = - V del " " . (11 



\n,n=l \ 




The G blocks denote green function matrixes at two sites and Go« = (Gjo)'- The short distance 
Green functions are given in Appendix, therefore to compute the above determinants at the scaling 
limit we have to know the expansion of Green function at far distances. This could be obtained by 
solving the deformed continuum Laplace equation: 
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where 7 = 0.577 ... is the Euler-Mascheroni constant, 
find the two point probability: 
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(12) 



Using all these together, we are ready to 



P(h < 1 - e, h r < 1 - e) = (P(l - 6)) 2 + f -^f L + ■■■ 



(13) 



where / is derived as an analytic functions in terms of e and 8. Again it is a very long expression 
and we do not bring it here. We have sketched the function / in polar coordinate for two different 
values of e = 0.1 and 0.2 in Figure 8. The result shows that the correlations are long ranged and 
have scaling property, though there is no rotational symmetry, and the amount of this anisotropy 
depends on e. 
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Figure 9: Three different configurations to modify the toppling matrix near a boundary 
4.3 Boundary effects on height 1 — e probability 

As boundary sites play important role in sandpile models, it would be interesting to consider the 
model in presence of different boundary conditions (BCs). We consider two types of boundary 
conditions, open and closed [26], according to toppling rules in our model. In the case of open BC, 
we set Ajj = 4 for boundary sites. This means that the amount of sand leaving the system when a 
boundary site topples is 1 — e or 1 + e depending on whether the site is on a horizontal or vertical 
edge. For closed BC we set An = 3 + e on horizontal edge and An = 3 — e for vertical edge so 
that there would be no dissipation on the boundary. It maybe hard to keep track of vertical and 
horizontal edges, an easier way is to consider e G (—1,1), so the transformation e — > — e would 
interchange vertical an horizontal edges. 

For simplicity, we consider that the model is defined on the upper half plane and the open/closed 
boundary is located at y = 1. Let's compute the probability of having an amount of sand less than 
1 — e at site % = (0, y) in the upper half plane. To this end we have to compute the open and closed 
Green functions for two sites i\ = (xi,yi) and % 2 = (£2,2/2)- This could be done using the image 
method: 

G op (x 1 ,y 1 ;x 2 ,y 2 ) = G(xi - x 2 ,yi - 2/2) - G{x\ - x 2 ,yi +2/2) (14) 
G d (xi, 2/15x2,2/2) = G{x\ -x 2 ,yi - 2/2) + G(xi - x 2 ,yi + 2/2 - 1) (15) 

The calculations of height 1 — e probability in the presence of the boundary are similar to that on 
the plane. There are three different configurations as shown in Figure 9, therefor^ 

P°P.<*(1 - €i y) = I (2 det (/ - B^G op f) + det (/ - B®G%f) + det (/ - B^G° h p /)) . (16) 

Using Eq. (I12p we are able to compute the determinants for distances far from boundary, then we 
find that the probability for height 1 — e at a distance y from boundary is: 

pop,d {1 _ e?y)=P{1 _ e) ±^ + (±\ (17) 

2/ \y J 

where the constant term P(l — e) is bulk probability for a site to have height at most 1 — e and is 
given by Eq. ([9]). The leading term that depends on the distance from boundary falls off as 1/y 2 

1 It turns out that the configurations B and C have the same contribution at large distances. 
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with coefficient 5(e) which is a complicated analytic function of e and is sketched in Figure 10. 
Plus sign corresponds to open boundary condition and minus sign to closed ones. It can be shown 
that this probability at distance x from a vertical open or closed boundary has a similar behavior. 

4.4 Field Theory Description 

As mentioned before, the number of different recurrent configurations of sandpile model (the parti- 
tion function) is given by determinant of the toppling matrix A which leads to the action of c = —2 
logarithmic conformal field theory (Eq. ©). We rewrite the action in x-y coordinates: 



S = j -9d z x 9 - 6d;e. (18) 

The lattice rotational symmetry is broken when we apply elliptic asymmetry in toppling rules. So 
it is expected the field theory describing our model not to be conformal invariant. The perturbed 
action turns out to be 

S~ J -ed 2 x 6 - 6d 2 y 6~ - e{ed 2 x 6 - 6d 2 y 6), (19) 

which is rotationally asymmetric. 

In this field theory, following the grassmanian method used in [II], it is possible to assign a 
scaling field to height variable being in [0, 1 — e). The probability of height 1 — e can be written as 

J Mi I dOj exp g O.X/lj) \ (exp (-9jB^) + exp (~^g%)) 

( e) " fdeifSjexpCEOi&ijOj) ' 

which resembles the expectation value of the field (1/2) fexp(— 9B^9) + exp(— OB^ffu . Following 
|14j . we can assign a field to this height variable in the following way: 

0(1 - e) = -((expt-^fl,) + expi-OiB^Sj))), (21) 

where ((. . .)) means that we have to contract all #'s and 9's except two. In the contraction we have 
to use the Green functions given in appendix. Doing all this and going to continuum limit we will 
arrive at 

0(1 - e) = - (f(e)d x 9dj + g(e)d y 9d y 9) (22) 
/(e) and g(e) are complicated functions with long terms that are plotted in Figure 11. 
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Figure 11: /(e) and g{e) are sketched as functions of e. 



As we see, setting e = we will arrive at the known field in ASM. As we increase e the function g 
gradually approaches to zero and is always positive. However the function / increases up to e ~ 0.2 
and then comes down becomes zero at e = 1. This is natural, because the limit e — ► 1 means you 
want to find the field associated with the probability of height h be in the interval [0, 1 — e] = 0, 
therefore the field should be zero. 

At any case it is very interesting to investigate this new field theory. In this theory, the scale 
symmetry is preserved and we expect criticality. However the rotational symmetry is broken and 
we do not have the full conformal invariance, though the system is still solvable from field theory 
point of view. Possibly on can investigate this model in the context of perturbed conformal field 
theory \27\ I28j . It can be argued that the action (|19p could become symmetric if you rescale 
horizontal and vertical directions with suitable (and different) scale parameters [29] . therefore this 
new theory should be identical with the previous one. This argument, if applied to the action is 
completely correct, and if one construct the field theory through defining the correlation functions 
the two theory will be identical. But if you ask questions about the field assigned to physical 
properties, such as the probability of height of a site being less that 1 — e, the situation becomes 
more complicated and it is not easy to find an easy way to relate the corresponding fields in the 
two theory. For example the field associated with the probability of height of a site being less that 
1 — e in ordinary CASM is (1 — e)(j)\ where <j)\ is the field of height one in ASM. However the same 
field in asymmetric case is given by equation (|22p and could not be derived from <j)\ via a simple 
rescaling. On the other hand it has been seen that the scaling exponent t s is not sensitive to the 
value of e. Putting all these together, it seems that the elliptical asymmetry do not change the 
universality class of the theory. 

It is also possible to restore the rotational symmetry statistically. If we assume the on-site 
asymmetries have quenched randomness; that is, on some sites e takes positive values and on some 
other it takes negative values, then there is no preferred direction statistically. The model should be 
formulated with more care, however it is very interesting to see if adding such quenched randomness 
will take the system to a new fix point or not. Work in this direction is in progress. 
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5 Appendix 

We collect some values of Green function G that we have used through the paper: 

1 
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